A First-Passage Kinetic IVIonte Carlo IVIethod for 
Reaction-Drift-Diffusion Processes 



Ava J. Mauro * Jon Karl Sigurdsson ^ Justin Shrake ^ Paul J. Atzberger § 

Samuel A. Isaacson ^ 



Abstract 

Stochastic reaction-diffusion models are now a popular tool for studying physical systems in which both the ex- 
plicit diffusion of molecules and noise in the chemical reaction process play important roles. The Smoluchowski 
diffusion-limited reaction model (SDLR) is one of several that have been used to study biological systems. Exact 
realizations of the underlying stochastic process described by the SDLR model can be generated by the recently 
proposed First-Passage Kinetic Monte Carlo (FPKMC) method. This exactness relies on sampling analytical solu- 
tions to one and two-body diffusion equations in simplified protective domains. 

In this work we extend the FPKMC to allow for drift arising from fixed, background potentials. As the correspond- 
ing Fokker-Planck equations that describe the motion of each molecule can no longer be solved analytically, we 
develop a hybrid method that discretizes the protective domains. Our discretization is chosen so that the dynamics 
of molecules within protective domains are approximated by continuous-time random walks. The empirical numeri- 
cal convergence and accuracy of our method is demonstrated in one dimension for both smooth and discontinuous 
potentials. As an application, we Investigate a model in which two molecules that can react undergo drift-dlffuslon 
along a polymer, but may also unbind from the polymer and diffuse in three dimensions. We study the interplay 
between unbinding rate, choice of potential, and polymer shape in the statistics of the reaction time and location. 

1 Introduction 

A fundamental challenge in cell biology is to understand how to predict and control the dynamics of cellular pro- 
cesses [2]. Stochasticity in the quantities and movements of molecules can have significant effects on the outcomes 
of cellular processes, particularly given the low copy number of many signaling and regulatory proteins and mRNAs 
present in a cell. For such species, the actual number and locations of molecules can provide a more accurate and 
useful description than the local concentration. The method we present in this paper allows for the explicit simu- 
lation of the stochastically varying numbers and locations of molecular species undergoing chemical reactions and 
drift-diffusion. 

Experimental studies and mathematical models have shown that stochasticity in the chemical reaction process 
plays a role in many cellular processes, for example gene expression [7, 14, 42], cell-fate decision making [47, 32], 
and signaling pathways in development [4]. Mathematical models of such processes frequently treat an individual 
cell as a single well-mixed volume or a small number of well-mixed compartments, such as the nucleus and the 
cytoplasm. However, the heterogeneous spatial distribution of chemical species, as well as interactions with internal 
membranes and organelles, often have significant effects on cellular processes. For example, after gene regula- 
tory proteins enter the nucleus through nuclear pores, the time required to find specific DNA binding sites can be 
significantly influenced by the spatial structure of DNA within the nucleus [27]. Similarly, the spatial distribution of 
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components of cellular signaling processes can play a decisive role in the successful propagation of signals from the 
cell membrane to the nucleus [36, 30]. 

A number of stochastic reaction-diffusion models have been introduced to understand the combined influence of 
noise due to the chemical reaction process and spatial diffusion [1,13, 48, 27]. These mathematical models resolve 
the explicit spatial movement of proteins and mRNAs within cells. There are three such mathematical models that 
have been commonly used: the spatially-continuous Smoluchowski diffusion-limited reaction model (SDLR) [46, 29]; 
what we call the Doi model [49, 10, 1 1]; and the lattice-based reaction-diffusion master equation model (ROME) [16, 
33]. 

These models often treat the movement of molecules as purely diffusive; however, drift can also play a significant 
role in the dynamics of cellular processes. Examples of sources of such drift include active transport, variations in 
chemical potential, material heterogeneities in the cytoplasm and nucleoplasm, and local interactions with cellular 
structures. The incorporation of drift has played a key role in developing models for molecular-motor based active 
transport [40, 5], movement of proteins on DNA [45], and protein movement subject to the influence of volume ex- 
clusion by chromatin [27]. It has been shown how to extend the ROME to incorporate drift due to potentials [27]. In 
the present work we consider a generalization of the SDLR model that allows for drift due to a fixed potential. For 
creating realizations of the stochastic process described by this model, we present a new numerical method com- 
bining elements of the First-Passage Kinetic Monte Carlo (FPKMC) method [12, 48] and the RDME-based methods 
of [52, 27]. 

1.1 SDLR and RDME 

In this subsection, we compare the SDLR and RDME models to provide motivation for our choice of the SDLR model. 
In the standard SDLR approach [46], the positions of molecules are modeled as point particles undergoing Brownian 
motion. The state of a chemical system is given by the collection of stochastic processes for the positions of each 
molecule of each chemical species at a given time. Bimolecular, or second-order, reactions between two molecules 
are modeled as occurring either instantaneously or a with fixed probability per unit time, when the molecules' sepa- 
ration reaches a specified reaction radius [29]. Bimolecular reactants are not allowed to approach closer than their 
reaction radius. Unimolecular, or first-order, reactions involving a single molecule represent internal processes, such 
as decay or splitting, and are assumed to occur with specified probabilities per unit time. (The Doi model is similar, but 
allows bimolecular reactants to approach arbitrarily close to each other. Bimolecular reactions then occur with a fixed 
probability per unit time when the separation between two molecules is less than the reaction-radius [49, 10, 11].) 
It is conceptually simple to extend the SDLR model to incorporate drift due to potentials. The underlying reaction 
process remains unchanged, but molecules move by a drift-diffusion process instead of pure Brownian motion. 

In contrast to the spatially-continuous SDLR model, space in the RDME is partitioned by a mesh into a collection 
of voxels. The diffusion of molecules is approximated by a continuous-time random walk on the mesh, with bimolec- 
ular reactions occurring with a fixed probability per unit time for molecules within the same voxel. Unimolecular 
reactions are modeled in the same manner as in the SDLR model. The state of a chemical system is then described 
by the collection of stochastic processes for the number of each chemical species within each voxel at a given time. 
Molecules are assumed to be well-mixed within each voxel (i.e. uniformly distributed). 

The standard RDME has one major drawback; in the continuum limit that the lattice spacing is taken to zero 
bimolecular reactions are lost in two or more dimensions [24, 26, 20]. This occurs because molecules are modeled 
by points, and as the lattice spacing approaches zero each voxel of the mesh shrinks to a point. Since two molecules 
must be in the same voxel to react, and in two or more dimensions two points can not find each other by diffusion, 
bimolecular reactions will never occur. Recently a modified convergent RDME (CRDME) that approximates the Doi 
model was proposed in [25]. While the CRDME gives a convergent lattice model, there are still several challenges 
to using such lattice models to study cellular processes in realistic geometries. Foremost, existing RDME methods 
approximate the domain geometry by either unstructured meshes [15] or Cartesian grid embedded boundary meth- 
ods [28]. In the former it can be difficult to construct meshes for which a continuous-time random walk approximation 
of diffusion is well-defined at all mesh voxels when in three dimensions [15]. The latter tends to lose accuracy in 
voxels cut by the boundary [28]. 

For these reasons, in this work we have chosen to focus on numerical methods for approximating the more 
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microscopic SDLR-based models. In particular, our approach does not constrain the molecules to remain on a fixed 
lattice as in the RDME. While our method is presented in one dimension, it can be adapted to avoid the issues 
described above for piecewise linear (planar) boundaries in two dimensions (three dimensions). 

Exact numerical realizations of the stochastic processes described by the RDME can be generated by the well- 
known Gillespie method [18], which is a Kinetic Monte Carlo (KMC) method [8]. More recently, the First-Passage 
Kinetic Monte Carlo (FPKMC) method was developed to generate exact realizations of the stochastic processes 
described by the SDLR model [38, 37, 12, 48]. It should be noted that there are a number of alternative numerical 
methods that have also been proposed for approximating the stochastic processes described by the SDLR or Doi 
models [3, 51 , 31 , 21]. Each of these methods is approximate, and ultimately based on a discretization in time. 

While it is conceptually simple to modify the SDLR model to include drift due to potentials, there are no standard 
methods for generating numerical realizations of the underlying stochastic processes described by the model. In 
this work we propose a modified version of the FPKMC method that allows for the inclusion of drift due to a fixed 
potential. While the method is no longer exact, the error introduced in the method is controlled. The method can be 
extended in a straightfonward manner to also include spatially varying diffusion coefficients. 

1 .2 FPKMC Method for the SDLR Model with Drift 

The central idea of the FPKMC method is to enclose one or two molecules within a 'protective domain' that isolates 
them from all other molecules of the system. A significant change in the state of the system occurs only when 
a molecule leaves its protective domain or a reaction occurs. In the cases discussed below, first-passage time 
distributions for such events can be determined from information about the underlying physical system. The first- 
passage time distributions can then be sampled to determine the time and type of state change that occurs next. 
This event-driven approach provides an especially efficient simulation algorithm by allowing each update of the 
algorithm to span the time interval between significant state changes (as opposed to proceeding over many small 
fixed time steps during which the change in relevant state variables is minor). 

In the rather special case of pure Brownian motion in simple domains (spheres or rectangular regions), the first- 
passage time distributions for a molecule to leave a protective domain [38, 37, 12, 48] or for two molecules to reach 
a threshold radius for reaction [48] can be computed analytically by solving the diffusion equation. The use of these 
expressions allows the generation of exact realizations of the stochastic process described by the SDLR model with 
the FPKMC. However, for many situations in cell biology, pure Brownian motion does not provide the most realistic 
description of the movement of molecules as a consequence of active transport, chemical gradients, interactions 
with cellular structures, etc. In such cases, significant drift terms are inherent to the particle dynamics and can 
be modeled as arising from a fixed potential. The method we develop extends the FPKMC to allow for such drift. 
Analytical expressions for the first-passage time distributions from protective domains are no longer possible with 
the addition of spatially varying drift. We therefore approximate the drift-diffusion process each molecule undergoes 
within a protective domain by continuous-time random walks on a discretized mesh. The "hopping rates" for these 
walks are derived by finite-difference discretizations of the Fokker-Planck equation. For this reason, our method 
can also be interpreted as a dynamic-lattice version of the RDME. Unlike the standard RDME, it has the benefit of 
converging to the SDLR model as the underlying lattice is reduced. 

The paper is organized as follows. Section 2 presents our approach for incorporating drift into the SDLR model. 
In Section 3 we give an overview of the implementation and steps of the FPKMC algorithm to generate realizations 
of the stochastic process described by the SDLR model. The general discussion in Sections 2 and 3 assumes that 
molecules move in M". The specific implementation we develop in Section 4 restricts molecules to intervals in M, but 
our technique can be extended to higher dimensions through the use of the "walk on rectangles" method [9]. Section 4 
presents our numerical method for incorporating drift into the FPKMC algorithm, and in Section 5 we demonstrate 
both the convergence and accuracy of the method in one dimension as the mesh spacing in the discretization is 
decreased. In particular, we apply our algorithm to the bimolecular reaction A + B ^ where the molecules 
of species A and B undergo drift-diffusion subject to various types of potential functions (smooth, discontinuous, 
and constant). Our results indicate that the method is approximately second-order accurate for smooth potentials 
and approximately first-order accurate for discontinuous potentials. Section 6 demonstrates several applications. 
In 6.1 we compare the effects of drift due to several potentials on reaction time and location statistics. We conclude 
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in 6.2 by investigating a simplified model of a coupled protein-polymer fiber system. The model takes into account 
polymer geometry and proteins that can undergo both one-dimensional drift-diffusion along the polymer and three- 
dimensional diffusion. 

2 Incorporating Drift into tlie SDLR Mo6e\ 

In this paper, we develop a First-Passage Kinetic Monte Carlo method for generating realizations of stochastic 
reaction-drift-diffusion processes, where the drift arises due to a potential. The stochastic reaction-drift-diffusion 
processes that our FPKMC method will simulate can be described by the modification of the SDLR model presented 
in this section. The general SDLR reaction-drift-diffusion model can be described by a system of partial integro- 
differential equations (PIDEs) for the probability densities of having a given number of molecules of each chemical 
species at a specified set of positions, similar to the stochastic reaction-diffusion PIDE models in [23] and [10, 11]. 
Alternatively, one can consider the collection of stochastic processes for the numbers of molecules and positions of 
each molecule of each chemical species in the system. Due to the high-dimensionality of the system of PIDEs for 
the probability densities, numerical methods for solving the SDLR model, including the FPKMC method, are typically 
based on Monte Carlo approaches that approximate the underlying stochastic processes. 

In the modified SDLR model that includes drift, molecules are modeled as points undergoing drift-diffusion pro- 
cesses. In a system with K chemical species, we label the k*^ chemical species by S'', k = 1, . . . ,K. We denote 
by M''{t) the stochastic process for the number of molecules of species S'' at time t. The position vector of the l^^ 
molecule of species S'^ at time t is given by the vector stochastic process Qf (t) e -R", Z = 1, . . . , M''. 

In the absence of any possible chemical reactions we assume the /"^ molecule of species S'' with position Qi(t) 
undergoes diffusion with diffusion coefficient D'' and experiences drift due to a potential V''(Qi{t)). In this case, 
Qi{t) satisfies the stochastic differential equation (SDE) 

dQtit) = ^ VF'=(Qf(i)) dt + V2D^dWf{t), (1) 

where A:b is Boltzmann's constant, T is absolute temperature, V denotes the gradient in the coordinates of Qi{t), 
and Wf (t) denotes the standard n-dimensional Wiener process which describes Brownian motion. Our method 
for generating realizations of this stochastic process will be presented in Section 4. The form of drift due to a 
potential in (1) is useful for modeling environmental interactions such as volume exclusion, attractive DNA forces felt 
by regulatory proteins, and effective potentials felt by molecular motors. It does not allow the possibility of interactions 
between the diffusing molecules of the chemical system. To incorporate such potentials into our method would be 
feasible for short-range pair interactions, but for simplicity we assume no potential interactions between molecules. 

When bimolecular and unimolecular reactions are added into the system, the molecules continue to move by 
the drift-diffusion process (1) under the additional constraint that any pair of bimolecular reactants are not allowed 
to approach closer than their corresponding reaction radius. Upon reaching this reaction radius we assume reac- 
tants either react immediately or with some specified probability per unit time. The former reaction mechanism is 
modeled through the use of an absorbing Dirichlet boundary condition, and the latter is modeled through the use of 
a partial absorption Robin boundary condition in the corresponding system of PIDEs [29]. Unimolecular reactions 
representing internal processes are modeled as occurring with exponentially distributed times based on a specified 
reaction-rate constant. We will assume that reaction products are placed at the locations specified in Table 1 . In 
the cases of two reaction products, for either unimolecular or bimolecular reactions, the angular orientation of the 
product separation vector about the center of mass is chosen randomly [3]. 

When using the Robin boundary condition-based bimolecular reaction mechanism, the separation distance will be 
chosen to be exactly the reaction radius. In the case that two products may react by the Dirichlet boundary condition- 
based mechanism, such a choice would lead to an immediate re-association reaction. This issue is generally handled 
through the introduction of an unbinding radius [3]. The need for unbinding radii can be avoided by using the Robin 
boundary condition-based reaction mechanism for all reversible bimolecular reactions. 
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Table 1 : Placement of reaction products. 





One Reaction Product 


Two Reaction Products 


Unimolecular 
Reaction 


The product is placed at the 
same location as the reactant. 


The products are placed a specified distance apart, with 
their center of mass at the location of the reactant. 


Bimolecular 
Reaction 


The product is placed at the 
center of mass of the two re- 
actants. 


The products are placed a specified distance apart, with 
their center of mass at the same location as the center 
of mass of the reactants. 



3 First-Passage Kinetic l\/lonte Carlo l\/letliods 

3.1 Overview of FPKMC Approaches 

First-Passage Kinetic Monte Carlo methods have been developed to generate exact realizations of the SDLR model 
in the absence of drift [38, 37, 12, 48]. These novel algorithms are based on the "walk on spheres" method for solving 
exit time problems in complicated geometries [35]. They rely on being able to derive exact analytical solutions of the 
diffusion equation in spheres and rectangular solids. In these FPKMC algorithms, a spherical or rectangular region 
called a 'protective domain' is drawn around every molecule in the system, with the collection of protective domains 
chosen to be disjoint. The first-passage time for each molecule, meaning the time when the molecule will first hit 
the boundary of its protective domain, can be sampled exactly using the corresponding analytical solution to the 
diffusion equation. The molecule that exits its protective domain first is updated to its exit position, and a new 
protective domain is defined. When two reactants are sufficiently close that a sphere can be drawn that contains only 
them, the corresponding two-body solution to the SDLR PIDEs can be used to exactly sample a candidate time and 
location for their reaction [48]. 

A central feature of each of these prior methods is the assumption that the molecular species undergo purely 
diffusive stochastic dynamics within their protective domains. When drift due to a potential is present in addition to 
diffusion, the probability densities for the locations of the molecules are no longer described by the diffusion equation, 
but rather by a Fokker-Planck equation. The general Fokker-Planck equation is given by 

= V • (p(x, t)VV{x) + Vp{x, t)) . (2) 

Here p{x,t) denotes the probability density the molecule is at x at time t and V{x) the strength of the potential 
at X. Our methods for using discretizations of several Fokker-Planck equations to sample times and locations of 
first-passage and reaction events are described in Section 4. Boundary conditions for the protective domains are 
also discussed in Section 4. 



3.2 Main Steps of the FPKMC Algorithm 

In this section, we describe the role of protective domains and the processing of events in our implementation of 
the FPKMC algorithm. We then list the main steps of the algorithm. Our implementation is based on the algorithm 
developed in [38, 37, 12], with some modifications. 

To apply the algorithm, every molecule in the system is placed in a protective domain. In one dimension the 
protective domains are intervals and in higher dimensions the protective domains are usually rectangular or spher- 
ical regions. In general, the boundaries of protective domains are absorbing. The boundary of a protective domain 
can contain a portion of the boundary of the overall spatial domain, in which case the protective domain boundary 
conditions will depend on the overall domain boundary conditions. We allow protective domains to contain either one 



6 




or two molecules. Protective domains containing only one molecule are referred to as 'single protective domains', 
and those containing two molecules are referred to as 'pair protective domains'. Molecules in separate protective do- 
mains behave independently. Each molecule undergoes drift-diffusion within its protective domain, and may undergo 
unimolecular reactions. Two molecules in the same protective domain may additionally participate in bimolecular 
reactions. To maintain independence when bimolecular reactants are in different protective domains, we require the 
boundaries of their protective domains to be separated by at least one reaction radius. For non-reacting molecular 
species we allow for overlap. 

Each event that may occur will have a type, time, and location. The two major event types are first passage from 
a protective domain and reaction. First passage from a protective domain occurs when a molecule first reaches an 
absorbing boundary of its protective domain. In the remainder we assume bimolecular reactions occur immediately 
when the separation between two reactants equals a specified reaction radius (a pure absorption reaction). Times 
and locations for first-passage events and for bimolecular reactants to first reach a separation of one reaction radius 
are sampled from drift-diffusion statistics, which are determined from the solution of the Fokker-Planck equation (2) 
(with appropriate boundary conditions to model the first-passage time problem for each protective domain). The time 
for a unimolecular reaction to occur is sampled from an exponential distribution with a specified reaction rate, and a 
corresponding location is sampled from the drift-diffusion statistics. 

To facilitate the discussion of these events, we use specific names for three times. The 'global time' will refer to 
when the most recent event has occurred, irrespective of its particular type or which molecules were involved. An 
'individual time' and a 'next event time' will be associated with each particular molecule; 'individual time' will refer to 
when the molecule was last updated, and 'next event time' will refer to the sampled time at which the molecule might 
next undergo an event. Individual times are less than or equal to the global time, and next event times are greater 
than the global time. 

Usually the individual time and location of a molecule will only be updated when the molecule undergoes a major 
event (first-passage or reaction). In this case, the time and location of the molecule will be updated to the time and 
location of the event. However, a molecule can also be updated to any specified time prior to its next event time, by 
moving the molecule to a location sampled at the specified time using the drift-diffusion statistics within its protective 
domain; this is referred to as a 'no-passage' update. 

The algorithm is carried out according to the following steps: 

1 . Protective domains are defined around each molecule or pair of molecules, as shown in Figure 1 b. 
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Figure 2: One simulation of the reaction A + B ^ 0, with one molecule each of A and B present initially and 
V{x) = 0. Both panels show the same run of the simulation. In the left panel, the vertical axis is the number of times 
the simulation cycled through steps 4 to 6 of the algorithm. In the right panel, the vertical axis is the time of the most 
recent event. 



2. The next individual event for each molecule or pair of molecules is determined by sampling an event type, time, 
and location. In Figure 1c, each next event time is labeled by a n. 

3. To determine global events, the individual events are sorted in a priority queue ordered from the shortest event 
time to the longest event time. For example, T5 denotes the shortest event time in Figure 1c. 

4. The next global event is determined from the priority queue using the next individual event with the shortest 
time. The global time and the individual time(s) of the participating molecule(s) are updated to the event time. 
In the case of a first-passage event for a molecule to leave its protective domain, the molecule's location is 
updated to the sampled first-passage location, as shown in Figure Id; if this molecule is in a pair protective 
domain with another molecule, the other molecule is no-passage updated to the new global time. In the case 
of a reaction event, the reaction products are placed at or about the reaction location, as specified in Table 1 . 

5. Molecules in protective domains that are close to the newly updated molecules are no-passage updated to the 
new global time. 

6. New protective domains are constructed only for those molecules that have undergone an update to reach the 
current global time, as shown in Figure If. New events are sampled for these updated molecules, and the 
event times are sorted into the priority queue. All other molecules and events remain unchanged. 

7. Steps 4 through 6 are then repeated. 

Note that Step 5 is used to keep the sizes of the protective domains from becoming too small [12], in which case the 
effective time steps used in the FPKMC method could become very short. 

We remark that Information about the state of any molecule in the system is available for any particular time in 
the simulation. For instance, if one would like to sample the locations of all molecules at a specified time, this can 
obtained by taking the state of the system at the largest global time before or equal to the specified time and then 
no-passage updating each molecule to the specified time. 

3.3 Protective Domain Clianges during One Simulation 

During simulations, updates are made to the protective domains sequentially as events occur changing the state of 
the system. To illustrate this process, we consider the simulation in one dimension of the reaction A + B ^ 
starting with one molecule each of A and B and using our FPKMC algorithm. One simulation is shown in Figure 2. 

In the left panel of Figure 2, the vertical axis is the number of times that the simulation cycled through steps 4 to 
6 of the algorithm; we call this number N. In the right panel, the vertical axis is the time of the most recent event. 
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At = 1 in this particular run of tlie simulation, molecule A is first-passage updated to the right endpoint of its 
initial protective domain. This location is close to the left endpoint of molecule B's protective domain, so molecule 
B is no-passage updated and new protective domains are defined around each molecule. From AT = 1 to iV = 3, 
molecule B is first-passage updated but does not come close to the protective domain of molecule A, so molecule A 
is not updated. At = 6, the distance between molecules A and B is less than a specified pair threshold, so they 
are placed in a pair protective domain. A\ N = 7, the distance between the molecules reaches the reaction radius 
and the reaction occurs. 



4 Implementation of FPKMC with Drift in One Dimension 

To simplify the discussion of our methods for introducing drift-diffusion into the FPKMC approach, we consider the 
case where the simulation domain is a one-dimensional interval. Also for simplicity, we assume that bimolecular 
reactions occur instantaneously when the reactants' separation reaches the reaction radius. We allow the simula- 
tion domain to have reflecting, absorbing, or periodic boundaries. Reflecting boundaries are modeled using zero 
Neumann boundary conditions, and absorbing boundaries are modeled using zero Dirichlet boundary conditions. 
Protective domains are proper subintervals of the overall domain. 



4.1 Propagation of Molecules within their Protective Domains: Incorporating both Drift 
and Diffusion 

For pure diffusion, the first-passage times, first-passage locations, and no-passage locations can all be determined 
from analytic solutions of the diffusion equation [38, 37, 12, 48]. In contrast, once drift is considered such analytic 
approaches are no longer possible in general. Instead one must consider probability densities that satisfy a Fokker- 
Planck equation such as (2), in which we restrict to the case where the drift arises from a spatially varying potential 
energy function V{x). To handle the FPKMC event statistics, we shall introduce approximations by treating the 
movement of each molecule within its protective domain as a discrete-space continuous-time Markov chain, more 
specifically a continuous-time random walk on discrete mesh points. Jump rates between neighboring mesh points 
are obtained using the spatial discretization of the Fokker-Planck equation from Wang, Peskin, and Elston [52]. We 
shall refer to this approximation as the 'WPE discretization.' 
The specific Fokker-Planck Equation we consider is 

dp{x,t) n^/an^ 



dt dx \ dx ' dx 

We remark that the factor {ksT)'^ is absorbed into the potential function V. The WPE discretization weights then 
determine the jump rates {i.e. probabilities per unit time) for molecules to hop from one mesh point to another. In 
particular, the jump rate for a molecule to hop from the mesh point Xi to a neighboring mesh point Xj, in the case of 
a uniform mesh of width h, is given by 

V{x^)-V{x,) 

h?eMV{x3)-V{xi)]-V ^' 

We shall extend this for non-uniform discretizations below. 

The discretization given by (4) converges at second-order for smooth potentials, and moreover can handle dis- 
continuous potentials [52]. The discretization also satisfies a discrete version of detailed balance at equilibrium, 
which helps reduce artificial drift due to numerical discretization errors [52]. Furthermore, this discretization is con- 
sistent with the standard second-order-accurate discretization of the Laplacian operator, in that aij converges to ^ 
as V{xj) -V{xi) approaches zero. 

This discretization approach can also be extended to higher dimensions. The jump rates in each coordinate are 
then given by (4). Moreover, it is not difficult to modify (4) to incorporate a spatially dependent diffusion coefficient 
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Figure 3: The top row shows the non-uniform sub-meshes that are defined when two molecules are first placed in a 
pair protective domain. The second (resp. third) row shows the uniform mesh that is chosen if molecule A (resp. B) 
hops first. The uniform mesh width, hp, is chosen to exactly divide tr. Note that hi + h2 = hp. The rates ooi and 
ao2 are given by Eqn. (5). The locations of the sub-mesh points are chosen so that the distances |yo-2;i|, \yo-X2\, 
\yi - xo\, and \y2 ~ xq\ are all exactly divisible by hp. After one molecule hops to a new point on its sub-mesh, the 
new distance between the two molecules will be one of these four distances. Then, a new mesh of uniform width hp 
can be defined so that both molecules lie exactly on mesh points. 



D{x). For example, in the case that D{x) is continuous, the constant D in (4) can be replaced by [D{xi) + 
Dix,)]/2 [52]. 

To make use of our discretization approach, a mesh is defined within each protective domain so that every 
molecule is located at a mesh point. Each molecule then undergoes a continuous-time random walk on the mesh. 
Exact sample paths of the molecules' random walks are generated using the Gillespie Method [18]. The resolution 
of this process can be adjusted depending on the size of the protective domain and on the desired trade-off between 
computational efficiency and accuracy. As such, our specific approach for choosing the mesh width and the locations 
of mesh points is described in more detail in the next subsection. 

Each protective domain either has absorbing Dirichlet boundaries (p = 0) at both endpoints, or an absorbing 
boundary at one endpoint and a reflecting Neumann boundary {dp/dx = 0) at the other. Any protective domain 
endpoints on the interior of the overall domain are absorbing, as are endpoints that coincide with an absorbing 
boundary of the overall domain. If one endpoint of a protective domain is located at a reflecting boundary of the 
overall domain that endpoint is made reflecting. 

For a newly constructed protective domain containing a single molecule, we determine the molecule's next event 
time by sampling an exact random-walk path for the molecule to hop on the mesh points until it reaches an absorbing 
boundary of the protective domain. For pair protective domains with two molecules, we perform random walks for 
each molecule until either: (i) one molecule reaches an absorbing boundary of the protective domain; or (ii) the 
distance between the two molecules is equal to the reaction radius tr. The mesh width for pair protective domains 
is always chosen to exactly divide r^ , so that the reaction occurs when the two molecules are exactly one reaction 
radius apart. The time that a molecule reaches an absorbing endpoint is the first-passage time, and the endpoint 
that the molecule reaches is the first-passage location. A no-passage location at any specified time before the next 
event time can be obtained by finding the last time in the sample path less than or equal to the specified time and 
taking the location of the molecule at that time. 

Non-uniform mesh cells are used when needed to conform to a boundary or to move molecules onto a uniform 
mesh where the mesh width exactly divides the reaction radius, as will be described in Subsection 4.2. Let xq be the 
initial location of a molecule on a non-uniform mesh, with xi and X2 denoting the locations of the neighboring mesh 
points in either direction. Note that we may have either xi < xq < X2 or X2 < xa < xi (see Fig. 3, top row). Let 
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Figure 4: Non-uniform mesli cell at an absorbing Dirichlet boundary of a pair protective domain (left panel), or 
a reflecting Neumann boundary of a pair or single protective domain (right panel). Arrows are shown only where 
the jump rates differ from those given by Eqn. (4) for uniform mesh cells. The non-uniform mesh cell at a Dirichlet 
boundary is chosen so that the boundary lies on a mesh point (i.e., at an endpoint of the mesh cell) and /idu- < hp. At 
a Dirichlet boundary the jump rates aoi and ao2 are given by the non-uniform rates (5) with hi = hp and /12 = h^ir- 
At a Neumann boundary the non-uniform mesh cell is chosen so that the boundary lies at its center and /incu/2 < hp 
or hs- in this case, there is no arrow from xq toward the boundary due to the Neumann boundary condition. The 
jump rate aoi, going away from the Neumann boundary, is given by (5) with hi = hp or hs, and h2 = /incu. 



hi = \xo - xi\ and /12 = No - X2\. The jump rates from xq to xi or from 2:0 to X2 are given by 

n 2g Vix,)-V{xo) ^ 

°' hjihi + h2)e^p[Vix,)-V{xo)]-l ' ■ 

The non-uniform discretization (5) is derived in Appendix A by modifying the WPE discretization of the Fokker-Planck 
equation [52]. In the case V = 0, our non-uniform discretization reduces to the non-uniform second-order-accurate 
discretization of the Laplacian operator at a Dirichlet boundary given in [17]. To our knowledge, Eqn. (5) gives a new 
discretization of the Fokker-Planck equation for non-uniform meshes. 



4.2 Choosing the Mesh within Protective Domains 

Single-molecule protective domains with absorbing boundaries are chosen to be symmetric about the location of the 
molecule. A maximum mesh width for all single protective domains, h^^^, is specified. Then, for each individual 
protective domain, the mesh width hg is chosen to be the largest value less than or equal to hf^"^ that exactly divides 
the distance from the molecule to either endpoint of the protective domain. The mesh points are then chosen so that 
the molecule and both endpoints of the protective domain lie exactly on mesh points. Having the endpoints lie on 
mesh points allows us to enforce the Dirichlet boundary conditions at the endpoints using the jumps rates (4) without 
modification. 

For single-molecule protective domains with one absorbing endpoint and one reflecting endpoint, the mesh width 
hg is chosen to be the largest value less than or equal to hf^^ that exactly divides the distance from the molecule 
to the absorbing endpoint. A mesh is defined so that the molecule and the absorbing endpoint lie exactly on mesh 
points. The mesh is uniform with the exception of one non-uniform cell used immediately adjacent to the Neumann 
boundary, as shown in Figure 4, right panel. 

In pair protective domains, the mesh width hp is a specified value chosen to exactly divide tr. Each time that two 
molecules are placed in a new pair protective domain, the initial distance between the molecules will not necessarily 
be divisible by hp. Rather than perturbing the molecules, non-uniform mesh cells are used to move one of the 
molecules, as shown in Figure 3, so that a uniform mesh of width hp can be defined with both molecules lying exactly 
on mesh points. Since this uniform mesh is chosen based of the locations of the two molecules, the endpoints of the 
protective domain may not conform with the mesh. In this case, one non-uniform mesh cell is used at each endpoint, 
which may be Dirichlet or Neumann (see Figure 4). 
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5 Convergence of the Method in One Dimension 

We perform convergence studies on simulations of tlie annihilation reaction A + B ^ 0, where the molecular 
species A and B undergo drift-diffusion subject to various potentials on the interval [0. 1]. Our results demonstrate 
both the convergence and accuracy of our method as the mesh widths, /i™'^'' and hp, in the discretization are 
decreased. We denote by M^(t) and M^{t) the number of molecules of A and B, respectively, at time t. In the 
first set of convergence studies only two molecules are simulated, M'^(O) = Af-^(O) = 1. Each simulation runs 
until the two molecules have reacted. A large number (10^) of simulations are performed in order to minimize the 
statistical error, so that the error due to the spatial discretization and the rate of convergence can be studied. In 
the next set of convergence studies, multiple molecules each of A and B are simulated, M^{0) = M^(0) = 10 
or M^(0) = M^(0) = 50, and each simulation runs until all the molecules have reacted. We will denote the i^^ 
molecule of species A by Ai, and the location of Ai at time t by Qf{t). Bj and Qf{t) are defined analogously. In 
the case M^{0) = M^{Q) = 1, we will drop the subscripts i and j. 

The convergence studies are performed using three different potential functions: (i) zero potential, V^eroix) = 
(which results in pure diffusion); (ii) a cosine potential with two energy wells, Kos(a;) = cos(47r3:); and (iii) a step 
potential with one step, 



The step potential is used to demonstrate that our FPKMC algorithm with the WPE discretization of the Fokker-Planck 

equation can handle discontinuous potentials. Note that adding a constant to any potential would not change the 
results, since the Fokker-Planck equation depends on the derivative of the potential but not the potential itself. In 
particular, any constant potential would produce the same results as V{x) = 0. 

In all the convergence studies, the length L of the overall domain is 1 unit, the boundaries of the overall domain 
are reflecting, and the diffusion coefficient £) is 1 unit^/sec for both A and B. The values used for the reaction radius 
tr will be specified in each subsection. We will use the notation ZY(a, b) for the uniform random distribution on the 
interval [a,b]. The initial locations Qf{0) and Qf{0) are drawn from U{a,b), where [a,b] ^ [0,L] will be specified 
in each subsection. If |Q^(0) - Qf*(0)| < tr for some i* and j* , then and Bj^ will react immediately, at t = 0. 

Since hp is always chosen to exactly divide tr, it necessarily follows that hp < tr. For single molecule protective 
domains we allow hf'^^ to be larger than tr. In the convergence studies, we set hf^^ = khp where fc > 1, and 
hold the ratio of h^'^'^ to hp constant as both are reduced to study convergence. As discussed in Section 4, the 
actual mesh widths hs used in single protective domains are almost always less than hf^"^, and non-uniform mesh 
widths are used in both single and pair protective domains. For these reasons, we keep track of the mean of the 
mesh widths that are actually used in each simulation. In calculating this mean, each mesh width is weighted by the 
number of times that it is actually used in a sample path. Then, for all simulations performed with fixed /i™^"" and hp, 
we calculate an overall mean mesh width by taking the arithmetic mean of the means for each simulation. 

In general, to describe a stochastic reaction-drift-diffusion system of many molecules by the probability density 
of having a given number of molecules at specified positions, a large coupled system of partial integro-differential 
equations is required [23]. In the special case of only two molecules, M^{0) = M^(0) = 1, with both molecules 
having the same diffusion coefficient D, the reaction-drift-diffusion system A + B ^ \n ^D can be described 
by a single 2D PDE: a Fokker-Planck equation (or a diffusion equation when V is constant). Let p{x,y,t) denote 
the probability density for finding molecule A at location x and molecule B at location y at time t. Then p evolves 
according the equations 




dp 
dt 



= DV ■ {pVV + Vp) 



onO < x,y < L, \x - y\ > 
on \x-y\= tr, 



p = 



(6) 



dp 

dr] 



= 



on X = 0, L ov y = 0, L. 



where i] = rj{x, y) denotes the outward pointing normal at the point [x, y). The 2D domain for Eqn. (6) is illustrated 
in Appendix B, Figure 17. When the initial locations of the two molecules in the FPKMC simulations are drawn from 
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U{0, L), the corresponding initial condition for the 2D Fol<l<er-Plancl< or diffusion equation is a constant, p{x, y, 0) = 
1/1/2. ^Q^Q^ jr, the following we define p on < a;, y < L by also defining p{x, y,t) = for \x - y\ < tr. 

Let T denote the random variable for the time at which the two molecules react. Using the solution p{x, y, t), we 
can calculate the survival probability, 

S{t) = Vv[T>t]= ( [ p{x,y,t)dxdy, (7) 
Jo Jo 

and the mean reaction time, 

/•OO /"OO 

E[r] = -/ tS'{t)dt^ S{t)dt. (8) 
Jo Jo 

Note that 1 - S{t) is the reaction time distribution function. The 2D Fol<l<er-Plancl< and diffusion equations (6) can 
both be solved numerically by finite difference methods, and the diffusion equation can be solved analytically using 
an eigenfunction expansion. These numerical and analytic solutions are discussed in Appendix B, and provide a 
baseline against which we compare the results of two-molecule FPKMC simulations in Subsection 5.1 . 

In what follows, we will use the term 'statistical error' to refer to the difference between the empirical value of a 
statistic (e.g., mean reaction time) estimated from the FPKIVIC simulations and the upper, or lower, bound of the 99% 
confidence interval for the statistic. By 'discretization error,' we will mean the difference between the empirical value 
from the FPKMC simulations and the actual value. In the two-molecule case, actual values are known exactly from 
the analytic solution, p(x, y, t), when V is constant, and are estimated from the numerical PDE solver described in 
Appendix B when V is not constant. 

Since we perform a large number of simulations (10^) at each mesh width in the two-molecule case, the statistical 
error is quite small, generally between 0.04% and 0.19% for various statistics. Although we perform fewer simulations 
(4 X 10^) when using multiple molecules each of species A and B, the statistical error is still reasonably small, 
generally between 0.4% and 1%. Our results show that as the mesh width is decreased, the discretization error 
rapidly decreases to below the statistical error. This demonstrates that our version of the FPKMC algorithm converges 
and accurately resolves the underlying reaction-drift-diffusion processes. 



5.1 Results of Two-Molecule Convergence Studies 

In this subsection, we consider the reaction A + B ^ ior a system with M"^(0) — M^(0) = 1. Here, the initial 
locations Q'^iO) and Q^(0) are drawn randomly from W(0, 1); the reaction radius is 0.02 units; and the pair 
threshold rpair is equal to 2rR (i.e, the molecules are placed in a pair protective domain when \Q^{t) - Q^{t)\ < 
''pair)- We study the convergence of the mean reaction time E [T] and the survival probability S{t) as the mesh 
widths h™"^"^ and hp are reduced. 

For V = 0, the errors in the FPKMC simulation results are calculated relative to the exact analytic solution, as 
determined from the eigenfunction expansion in Appendix B. For Vcos and 14tep. the errors are relative to the numer- 
ical solution from the PDE solver described in Appendix B. The mean reaction time calculated from the numerical 
PDE solution was resolved to an absolute error tolerance of 10~^ for Vcos and 10~^ for Vstep- 

Let Eempl?"] be the empirical mean reaction time calculated from the FPKMC simulations and define Eupp[T] to 
be the upper bound of the 99% confidence interval for the empirical mean reaction time. We denote by Eactl?"] the 
exact analytic mean reaction time in the case that y = 0, or the mean reaction time determined from the numerical 
PDE solution in the y ^ cases. We calculate the relative error by 

|Eact[r]-Eeinp[r]| E„pp [f] - Eemp [T] 

Eact[T] Eact[T] • ^ ' 

For each mesh width, the empirical survival probability and the associated 99% confidence bounds are calculated 
using the MATLAB function 'ecdf '. Let S'emp(i) be the empirical survival probability Supp{t) the upper bound of the 
99% confidence interval for the empirical survival probability, and Sact {t) the analytic or numerical survival probability. 
The discrete L^, L^, and L°° norms of a function u{t) are given by 

II^WIUi =IIw(^OAii, \\uit)\\L2 ^ [J2u{t^f^t^}\ | I |l~ = max (10) 
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Figure 5: Convergence of mean reaction time, E [T], for the two-molecule A + B ^ reaction, as the mesh 
width is decreased. Note the expanded scales of the vertical axes. Each FPKMC data point shows Eo„ip[T] from 
lO*" simulations. Error bars give 99% confidence intervals. In the = case, we compare the FPKMC simulation 
results to the analytical solution, as opposed to the Crank-Nicolson PDE solution for V ^ 0. The Crank-Nicolson 
results are shown here in the F = case to demonstrate the accuracy of the PDE solver discussed in Appendix B. 



The relative errors in Scmp{t) we report in each norm are then given by 

||'5'act(^) ~ ■5'cmp(^)|| I ||'5'upp(t) — Scmp{t)\\ . , 

ll^actWII ll^actWII ' 

where the norms are evaluated on the interval t e [0,ts]. Here is the time at which Sact{ts) = 10^*^, and 
Ati — ti+i — ti. The time points, tj, used in evaluating the norms correspond to those at which the numerical PDE 
solutions were calculated (see Appendix B for more information). In the = case, the analytic expression for the 
survival probability was evaluated at those same U's. In all cases, the empirical FPKMC survival probabilities were 
linearly interpolated to obtain values at every t^. 

In each of the Figures 5 through 7, the left panel shows results for V = 0, the center panel shows results for 
V"cos. and the right panel shows results for V^stop- The insets in Figure 5 show the respective potentials. In Figure 5, 
Eemp[r] is plotted against /if'' as the mesh width is varied. Figures 6 and 7 shows the relative errors in Eempf?"] (9) 
and in Scmp{t) (11), respectively, plotted against the mean mesh width. As the mesh width in the FPKMC simulations 
is decreased, EcmpiT"] approaches Eactl?"] and Scmp{t) approaches 5'act(0 for all three potentials. 

For = 0, the relative errors in Eomp[r] and S'omp(i) for each mesh width are all less than 1.13% and decrease 
to less than 0.21% for the finer mesh widths (Figs. 5-7, left). For V^os, the relative errors are all less than 6.57% 
and rapidly decrease to less than 0.22% (Figs. 5-7, center). Furthermore, the discretization errors decrease to 
less than the corresponding statistical errors, for both V = and Vcos- Note that the statistical errors are very small 
since we performed lO'' simulations at each mesh width. The results in the y = case show that the reaction and 
diffusion processes are very well resolved by our method. The results in the Vcos case demonstrate that our method 
accurately resolves drift due to a smooth potential, in addition to reaction and diffusion. For Vstep, the relative errors 
are less than 4.37% for all mesh widths and decrease to less than 0.38% (Figs. 5-7, right), which demonstrates 
that our method can handle discontinuous potentials. 

We estimate the rate of convergence of our method to be approximately second-order for V"cos and approximately 
fist-order for Knop- This is consistent with the convergence rates of the WPE discretization of the Fokker-Planck 
equation for smooth and discontinuous potentials. In the V" = case, it is difficult to draw a conclusion about the 
rate of convergence since the discretization errors are small relative to the statistical errors; however, for the same 
reason, we can conclude that our method is very accurate in this case. 
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Figure 6: Relative errors in the empirical mean reaction time, Ecinp[r], for the two-molecule A + _B — > reaction. 
Note that the vertical axes have different scales in each panel. Each data point is based on 10'^ simulations with 
fixed hf^"^ and hp. The statistical error bars are determined by Eqn. (9) using the 99% confidence intervals for the 
empirical mean reaction times. 




10-^ 10-1 10-2 10-1 10-2 10-1 
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Figure 7: Relative errors in the empirical survival probability, S'cmp(^), for the two-molecule A + B ^ reaction. 
Each empirical survival probability function is based on 10^ simulations with fixed hf^'^ = 8hp. The statistical error 
bars are determined by Eqn. (1 1 ) using the 99% confidence bounds for the empirical survival probabilities. 
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Figure 8: Mean time for all molecules to react via the reaction A + B ^ 0. Error bars show 99% confidence 
intervals, based on 4 x lO"* simulations per data point. Left panel: j-r = 0.02 units. Right panel: tr = 0.001 units. 



5.2 Results of Multiple-Molecule Convergence Studies 

In the convergence studies for multiple molecules undergoing the reaction A + B — > 0, we start the simulations 
with either 20 molecules (M^(0) = M^{0) = 10) or 100 molecules (Af^(O) = A/^(0) = 50). In the 20-molecule 
simulations: tr = 0.02 units; rpair = 2rR; Qf (0) - Z^(0. 1,0.4) and Qf (0) - W(0.6,0.9), < i, j < 10. In 
the 100-molecule simulations: = 0.001 units; rpair = 4rR; (0) and Qf (0) - U{0, 1), < i, j < 50. Two 
molecules Ai* and Bj^, are placed in a pair protective domain if they are closer to each other than to any other 
molecules of the opposite type, i.e. if 

\Qtit) - QfM\ = min \Qtit) - QfM = min \Qtit) - gf 
and if the distance between them satisfies 

\Qt{t)~QfAt)\<rp.ir 

and 

\QtAt)-QfAt)\<rR + rnm{mm{\Qt{t)-Qf{t)\, min \Qf, {t) - Qf {t)\). (12) 

The last condition (12) was added to prevent the protective domain of a molecule, say Ai, from becoming too small 
when Ai is close to another molecule of the same type, say Ai^, where would otherwise have been placed in a 
pair with B^*. 

Each simulation runs until all of the molecules have reacted. Let r„ denote the random variable for the time 
at which the n^^ reaction occurs. Since we are working with the irreversible reaction A + S 0, we have that 
M^(i) = M'^(O) — n for t e [T„,r„+i), and similarly for M^{t). Figure 8 shows the empirical mean time for all 
molecules to react; this is Eemp[rio] when M^(0) ^ M^{Q) = 10 and Eemp[?5o] when M^(0) = M^(0) = 50. 
With these values of M^{Q) and A/^(0), the SDLR model described in Section 2 will correspond to a large coupled 
system of partial integro-differential equations for the probability densities of having a specified number of molecules 
at specified locations. It is no longer feasible to solve these equations numerically to obtain high-accuracy solutions 
for assessing the empirical convergence of our FPKMC method. As such, we now estimate the accuracy of reaction 
time statistics by comparing FPKMC results at coarser mesh widths to results obtained with the finest mesh width. 

When V = cos(47ra;) (Fig. 8, top panels), we see convergence as the mesh width is decreased. The percent 
difference between Ecmp[Tio] for the coarsest mesh width compared to the finest mesh width is approximately 5.4%; 
this percent difference for Eomp[75o] is approximately 7%. These differences are comparable in size to the explicitly 
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Figure 9: Mean number of molecules of A remaining at time t, Ecinp[Af'^(t)] (left panel), and probability that at 
least one molecule of A remains at time t, Pr [M'^(i) > l] (center and right panels). V = cos(47ra;). Each graph is 
based on 4 x 10"* simulations. Eoinp[M^(i)] is not plotted in the case Af^(O) = M^{0) = 50, because the results 
for the different mesh widths are essentially indistinguishable until E^j^plM^it)] < 4. Left and center: tr = 0.02 
units, = &hp. Right: tr = 0.001 units, hf'''-' = 200hp. 

calculated discretization error of approximately 6.4% at the coarsest mesh width in the two-molecule case (previous 
subsection). 

When V = (Fig. 8, bottom panels), the confidence intervals for all mesh widths overlap, indicating that the 
discretization error is less than the statistical error even for the coarsest mesh width. The statistical errors when 
= are between 0.49% and 0.81%. If the unknown discretization error here is comparable in size to the known 
discretization error in the two-molecule case, which was approximately 0.87% at the coarsest mesh width, that would 
provide an explanation for why can not observe convergence when V = 0. 

Ecnip[M^{t)] is the mean number of A molecules remaining at time t, and Pr [M^{t) > l] is the probability that 
at least one molecule of A remains at time t. Figure 9 shows convergence of Ecmp[M'^{t)] and Pr [M^{t) > l] 
as the mesh width is decreased, in the case V = cos(47ra:). Results are not plotted for the y = case, because 
the confidence bounds for different mesh widths overlap; this indicates that the results are resolved to within the 
statistical error, even for coarser mesh widths. 

6 Applications 

6.1 Comparison of Potentials 

To demonstrate the contrasting effects that can be produced by different drifts, we consider reactions A + B ^ 
where the molecules diffuse within various potential energy landscapes. We consider the following three cases: (i) 
zero potential, V{x) = 0; (ii) a one-well potential, V{x) = cos(27ra:); and (iii) a two-well potential, V{x) = cos(47ra;). 
We use a domain of length L = 1 with reflecting boundaries, and diffusion coefficient D = I unit^/sec for both 
molecular species A and B. 

In the absence of reactions, the equilibrium probability density for a molecule to be at location x is given by the 
Boltzmann distribution 

We compare the Boltzmann distributions for each of the three potentials to the reaction locations from the A+B 
FPKMC simulations in the particular case of two molecules, A/'^(0) = A/^(0) — 1. The results are shown in 
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Figure 10: Reaction locations from A + B ^ FPKMC simulations with A'f^(O) = M^(0) = 1, and Boltzmann 
distribution for each potential. The reaction location is {Q'^it) + Q^{t))/2, where t is the time of the reaction. 
Each graph of reaction locations is based on 10 simulations, tr = 0.02. hf^^ = hp ^ 7'r/8. Q^{0) and 
Q^{Q) ^ U{0, 1). The plotted densities were determined by binning the reaction locations into 100 bins. 




Figure 1 1 : Left and center panels: Empirical survival probability S'cmp(i) for the two-molecule A + B ^ reaction. 
Both panels show the same survival probabilities, however the axes have different scales. The dashed lines show 
99% confidence bounds based on 10"" simulations, j-r = 0.02. hf^"^ ^ hp ^ tr/S. Right panel: Empirical mean 
number of molecules of A remaining at time t, for the reaction A + B ^ with M'^(O) — M^{0) = 10. The 
dashed lines show 99% confidence bounds based on 4 x 10"* simulations, tr = 0.02. hp = tr/IG. /i'"""" = 8hp. 



Figure 10. The potentials serve to spatially "confine" molecules, in the sense that molecules are most likely to be 
found in locations where the potential energy is lowest. Consequently, the reactions are most likely to occur in such 
low energy locations. 

As would be expected, the mean reaction time for the two-molecule A + B ^ reaction is faster with the 
one-well potential 0.03481 sec) than with no potential 0.06483 sec), while slower with the two-well potential 

0.09887 sec). Figure 1 1 (left and center panels) compares the survival probabilities S'cmp(i) for the three different 
potentials. The semi-log graphs of Scmp{t) for all three potentials appear linear except at short times, indicating that 
the reaction time distributions could be described by an exponential distributions for larger times (Fig. 1 1 , left panel). 
However, at short times (Fig. 11, center panel), the graphs are not linear and behave differently from each other. 
Initially, reactions occur more quickly with the two-well potential than with either the one-well potential or no potential. 
However, after 5cmp(^) has decreased to about 50%, reactions occur more slowly with the two-well potential than 
with the other potentials. As would be expected, if the initial locations of the two molecules are in the same well of a 
potential, then they tend to react more quickly; whereas, if the two molecules start in different energy wells, then the 
time until they react tends to be longer. 

Figure 11 (right panel) shows Ecinp[M^{t)], the mean number of molecules of A remaining at time t, when 
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Figure 12: Biopolymer diffusion-excursion search and protein-protein interactions. The geometry of an individual 
biopolymer or network is expected to be influential in the search kinetics of proteins which diffuse in one dimension 
along the biopolymers, but may also detach to undergo excursions in three dimensions. 

M^(0) M^{0) ^ 10, Qf{0) - U{0.1,0A) and Qf (0) - W(0.6,0.9). In this case, Ec^p[M^{t)] could be 
described by exponential distributions for all three potentials. In the case of A/^(0) = Af^(O) = 50 with Qf{0) and 
Qf{Q) ^ U{0, 1), Ecmp[M^{t)] is not plotted but behaves very similarly to S'omp(i) in the two-molecule case with 
Q^(0) and Q^{0) - U{0, 1) (Fig. 1 1 , left and center panels). 

6.2 Role of Biopolymer Geometry in Protein Diffusive Search 

Many proteins diffuse effectively in one dimension by adhering loosely to biopolymer filaments, such as actin, 
microtubules, or DNA [2]. In the case of DNA regulatory proteins, single molecule experiments and theoretical 
work indicate that observed rapid kinetic rates are achieved by a search process involving a combination of one- 
dimensional diffusion of the protein sliding along the DNA in conjunction with three-dimensional diffusive excur- 
sions [53, 43, 44, 19, 45, 34, 22, 6]. The geometry of the individual biopolymer or network of biopolymers is expected 
to play an influential role in such search kinetics. We introduce a model to investigate this phenomena in the diffu- 
sion of proteins that interact with biopolymers having basic knotted closed-loop configurations, see Figure 12. In our 
model we also allow for position-dependent forces acting on the proteins, such as might occur from heterogeneities 
in DNA sequence or protein binding sites. We perform simulations of two proteins undergoing a search process until 
they encounter each other on a biopolymer. The process involves one-dimensional drift-diffusion along the polymer 
in combination with absorption/desorption events associated with excursions from the polymer. For example, this 
process could model the formation of a regulatory complex at a non-specific DNA binding site. Movement along the 
polymer is simulated using the drift-diffusion FPKMC method we introduced in the preceding sections. The effects 
of the biopolymer geometry are taken into account through the absorption/desorption statistics computed by solving 
a diffusion equation in three dimensions for each specified biopolymer configuration. We remark that the model and 
methods we introduce are readily extended to more complex geometries, polymer networks, and filament bundles. 

6.2.1 Model of the Biopolymer Drift-Diffusion Process and Thiree-Dimensional Excursions 

The biopolymer is represented geometrically as a one-dimensional filament embedded within a three-dimensional 
space. The heterogeneity along the biopolymer strand, for example changes in DNA sequence, density of bound 
proteins, or other chemical factors are taken into account through the potential energy term in the drift-diffusion 
process. To account for diffusive excursions in three dimensions, we introduce desorption and absorption events 
for the protein from the biopolymer. Let Aoff be the first-order rate at which a protein desorbs from the biopolymer. 
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We assume the protein must undergo a conformational cliange before regaining affinity to bind to the biopolymer 
(otherwise it would instantly re-absorb at the desorption position). This conformational change is considered as a 
first-order reaction with the rate ^la- To model this process, let p(s2,t2|si,<i) denote the probability density for a 
protein that desorbs at location si at time ti to re-absorb at location 32 at time t2- 

This probability density can be calculated by considering the two stages of the protein desorption and re- 
absorption to the biopolymer. In the first stage, the protein desorbs at location si along the biopolymer and performs 
three-dimensional diffusion until the occurrence of the conformational change at a time Tq after desorption. At this 

time the probability density the protein is at the location X is given by ps(X \si,Ta) = (2^Df)-V2 "^^P 
where Z(si) denotes the position on the biopolymer from which the protein desorbed and D the three-dimensional 
diffusion coefficient of the protein. In the second stage, the protein has affinity for the biopolymer and diffusively 
re-absorbs to the biopolymer at the location S2 at time t2. We denote this probability density bypa(s2,t2|X,ti -j-Ta). 
This gives the density J3(s2,t2|si,ii) = // Pa(s2, i2|x, + f)ps(x|si, i) /Xae"'^"* dx. 

For simplicity we shall assume the biopolymer geometric conformation Is fixed, resulting In a desorption and 
absorption process that is stationary. This gives Pa{s2,t2\si,ti + Tq) = Pa{s2,t2 - ti - Ta\si,0). In practice, 
we tabulate pa(s2,t|X,0) at select locations, X, as a one-time off-line calculation using a numerical diffusion 
equation solver. This provides a very efficient method to simulate the excursions. To obtain Tq, we first generate 
a random exponential time at which the affinity conformational change of the protein occurs. We then sample the 
protein location, X = x, using the density Ps(X|si,Ta). From our pre-constructed table, we can then sample 
Pa{s2,t2 - ti - Ta\x) to find a time of re-absorptlon, t2, and the protein association location, S2- 

This provides an efficient method to simulate the diffusive excursions of proteins and their re-absorption to the 
biopolymer. The modeled search process of the two proteins proceeds by drift-diffusion along the biopolymer until 
either a desorption event occurs or the proteins encounter each other. The drift-diffusion along the biopolymer Is 
handled by our one-dimensional FPKMC algorithm with periodic boundary conditions. Upon a desorption event at 
location si, the protein is repositioned on the biopolymer at location S2 at the re-attachment time, t2, according 
to the density p(s2,t2|si,ti)- This repositioning process models the three-dimensional diffusive excursion of the 
protein until re-absorbing to biopolymer. The distribution p takes into account the role played by the geometry of 
the biopolymer conformation. An important feature of our model is the protein search process is reduced solely to a 
one-dimensional description. We remark that reflecting Neumann boundary conditions In our FPKMC method could 
be used to model obstructions on the biopolymer. Either a potential, or an absorbing DIrlchlet boundary condition, 
could also be used to model irreversible binding sites on the polymer. 



(X-Z(si))'' 



6.2.2 Numerical Methods for Excursion-Time Probability Distribution and Reattachment Locations 

We numerically solve the three-dimensional diffusion equation with diffuslvlty D for the position of a detached 
molecule at a given time. Denote by p(x, t) the corresponding probability density that solves the diffusion equation. 
Reattachment sites on a biopolymer are modeled using sink terms in the equation. We let xq label the detachment 
position, and s the parameterization variable along the biopolymer. From the probability mass absorbed at the loca- 
tions of the sinks, we can obtain the probability densities p{s, t\xo, 0). In the diffusion equation we use a sink term 
of the form 

g{y,t) = -(^J \{s)q{y-x{s))ds^ p{y,t), (13) 

where A is the intensity of the sink and q is an averaging kernel function. The rate of the reattachment flux is given 
by 

f (^'*) = / 7W^»^(^'*^^ XisHy-xismy,t),y. (14) 

In the discretized form this becomes 



aiVm^ = - X] '^kQiVm - x{Sk))ASkp{ym, t) (1 5) 

k 
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Table 2: Diffusion-excursion search parameters. 



3D Parameters 


Value 


1 D Parameters 


Value 


3D diffusion coefficient 


2.183823 Atm^sec-i 


1 D diffusion coefficient 


0.01/im^sec~^ 


Mo 


5.1728 X lO^ns-i 


Aoff 


0.02 to 200 sec"^ 


-^absorb 


5.00 X lO-^ns-i 


r-R 


20 nm 


Ax 


12.5nm 


/^max 


10 nm 


At 


4292.9ns 


hp 


5 nm 


-^biopolymer 


100 


^pair 


4 tr = 80 nm 
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and the rate of the flux becomes 

dQkjt) 
dt 



X^Afe5(y„-Xfe)Asfep(y„,t)Ay„. (16) 



The particular choice that we make for the kernel q is the four-point Peskin-5 kernel [39]. For x = {x, y, z) this is 
given by q{x) = (j)p{x/ Ax)(j)p{y/ Ax)(j)p{z/ Ax), where 

^ f 3 - 2u -t- Vl + 4u - Av? 0< «<1, 
(l}p{u) = j7.{ 5 - 2u - \/-7 -I- 12u - 4u2 i < u<2, 
i 2< «, 

with (l).p{-u) = (j)p{u). This particular function is chosen to reduce numerical error induced by the off-lattice shifts 
of the polymer adsorption locations relative to the underlying discretization lattice used in the three-dimensional 
diffusion solver. We take the absorption parameter A to have two stages. Namely, a "non-sticky" state A = and a 
"sticky" state A(.s) = Afe = Aabsorb- Note the latter is taken to be uniform along the biopolymer. As discussed in the 
last subsection, the protein becomes sticky after undergoing a conformational change, the time for which is modeled 
as an exponentially distributed random variable. 

We consider three distinct geometric configurations for the biopolymer, each with the same arc-length but a 
different topology Biopolymers represented by a circular unknotted loop, a trefoil knot, and a figure-eight knot are 
studied. We use parameterizations given by X{s) = c* {x{s),y{s),z{s)) for < s < 1, where c is chosen 
so that the arc-length of each of the knots is Lbiopoiymor = lOOOnm. The x, y, and z parameterizations used 
for the three polymer configurations are: unknotted loop, x{s) = cos(27rs), y{s) = sin(27rs), z{s) = 0; trefoil 
knot, x{s) = {2 + cos(67rs)) cos(47rs), y{s) = {2 + cos(67rs)) sin(47rs), z{s) = siii(67rs); and figure-eight knot, 
x{.s) = (2 + cos(47r.s)) cos(67r.s), ?y(.s) = (2 + cos(47r.s)) sin(67r.s), z{s) = sin(87rs). For each configuration we use 
the diffusion equation solver to tabulate the joint reattachment time and location distribution, see Figure 13. 



6.2.3 Simulation Results: Diffusion-Excursion Search with Different Biopolymer Geometries 

We assume one molecule each of protein species A and B are initially present on the biopolymer. The initial 

positions are drawn from a uniform distribution over the length of the biopolymer. The two molecules undergo drift- 
diffusion along the biopolymer and diffusion in the three-dimensional space by the method developed in the preceding 
subsections. Detachment times are sampled from an exponential distribution with rate Aoff. A reaction between the 
two molecules can only occur when both are on the biopolymer. 

For the parameters we take a = 10 nm to be the radius of the protein. The three-dimensional diffusion coefficient 
D = ksT/'j, where 7 = Gwrja is the Stokes drag, T = 298.15 is room temperature, /cb is Boltzmann constant, and 
the dynamic fluid viscosity, 77, is chosen to be 1 times the viscosity of water. The overall spatial domain is a box 
with sides of length 500nm. The parameter for the exponential distribution that represents when the protein becomes 



21 




Figure 13: Reattachment probability distributions for the different topologies, (a-c) Probability density of absorbing 
at a given location Z(s) when the protein desorbed at location Z(0.5). Ps^is) = p{s,t\si = 0.5,0) dt. (d-f) 
Conditional probability densities for a protein that desorbed at location Z(0.5) to absorb at a point S2. For select 
locations these are indicated in the first row by a colored dot with the same color convention holding for the plots 
in the second row. Psi.sa (t) = p{s2, t\si = 0.5, 0). (g-i) Biopolymer geometries from left to right: a circular loop, 
trefoil knot, and figure-eight knot. 




Figure 14: Reaction location densities \ox A + B ^ with M'^(O) = M^(0) — 1 and V{s) = 0. The reaction 
location is the midpoint between the locations of A and B at the time of the reaction. Each graph is based on 10^ 
simulations. The plotted densities were determined by binning the reaction locations into 50 bins. 
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Trefoil, V{s) = Qcos(67rs) "l/(s) = 1.5cos(67rs) 




Reaction Location (/im) Reaction Location (/im) 

Figure 15: Reaction location densities for A + S ^ with Af'^(O) = M^{0) = 1 and V{s) = acos (Btts). Each 
graph is based on 10^ simulations. The plotted densities were determined by binning the reaction locations into 50 
bins. 



sticky is denoted by fia = SL*/ y \ -^b^pp^y"" j . por convenience, the absorption strength, Aabsorb = 0.005ns ^, 
is chosen to be larger than the timescale required for numerical stability in the diffusion equation solver. The other 
parameters are specified in Table 2. 

For each of the two knotted biopolymer conformations, the reaction location distribution was investigated as the 
biopolymer detachment rate, \os, was varied, see Figure 14. These simulations assumed V{s) = so that there was 
no drift due to a potential. As Aoff increases, the effects of the trefoil and figure-eight biopolymer geometries become 
more pronounced. For the unknotted loop the reaction locations are approximately uniform over the biopolymer 
regardless of the choice of Aos. 

A natural question is how diffusion in a potential energy landscape may enhance or counteract the effects of 
biopolymer geometry on the reaction locations. From Subsection 6.1, we expect the density of reaction locations 
to decrease in areas where the potential is large and to increase in areas where the potential is small. To test this 
idea, we investigated the trefoil knot conformation with potentials of the form V{s) = acos (67rs) for several values 
of a e [-6,6]. 

The resulting reaction location densities as a is varied are shown in Figure 15 (left panel). As expected, the 
potential enhances the effects of the trefoil knot geometry when a < 0, and counteracts the effects when a > 0. 
When a = 1.5, the effects of the potential and the trefoil geometry essentially cancel each other out. For comparison, 
we also ran simulations with a = 1.5 for each of the other two biopolymer geometries, see Figure 1 5 (right panel). 

As the detachment rate Aoff is varied, the mean reaction times attain a minimum at approximately Aoff = 20, 
see Figure 16 (left panel). The mean reaction times are slightly faster with the trefoil and figure-eight biopolymer 
conformations than with the circle conformation. As Aoff is increased the mean total time that both molecules are on 
the biopolymer decreases monotonically, as shown in Figure 16 (center panel). 

When the potential V{s) — acos (67rs) is used with the trefoil biopolymer conformation, the mean reaction time 
is fastest when the amplitude |a| is large, see Figure 16 (right panel). As discussed in Subsection 6.1, molecules 
that are in the same energy well of a potential tend to react more quickly. However, if the proteins can only move on 
the biopolymer, potentials with large amplitudes as used here would represent large energy barriers, greatly slowing 
the time for the two molecules to find each other. For example, with a = -3 and Aoff = 0, the mean reaction 
time is approximately 80 ± 9 seconds, based on 1000 simulations (compare to Figure 16, right panel). The three- 
dimensional diffusion excursion provides an alternative path for proteins to circumvent the energy barriers present on 
the biopolymer. This could be an important mechanism in protein-protein interactions associated with biopolymers. 

Overall, these results demonstrate the importance of incorporating drift-diffusion and other spatial heterogeneities,! 
such as biopolymer geometry, when investigating biological chemical processes. They illustrate the power of the pre- 
sented FPKMC method in capturing such effects. Our approach and methods can be readily extended to more 
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V{s) = V{s) = 0, Circle V{s) = a cos(67rs), \os = 20 sec" 




Aoff(sec 1) Aofi(sec ^) a 

Figure 16: Mean reaction time and mean time tliat botli molecules are on the biopolymer, for A + i? — ^ witli 
M^{0) = A/^(0) = 1. Each data point is based on 10^ simulations. Error bars indicate 99% confidence inter- 
vals. The center panel shows results for the circle conformation only, but the results for the trefoil and figure-eight 
conformations are similar. 



complex chemical kinetics, biopolymer geometries, or filament networks. 

7 Conclusion 

We have presented a new First-Passage Kinetic Monte Carlo method (FPKMC) that is capable of incorporating the 
important roles played by drift and spatial heterogeneities into the stochastic dynamics of chemically reacting molec- 
ular species. We have shown that our numerical method is convergent for smooth potentials with approximately 
second-order accuracy and for discontinuous potentials with approximately first-order accuracy. Unlike the stan- 
dard RDME, our method retains bimolecular reactions as the protective domain lattice spacings are taken to zero, 
and converges to the underlying SDLR model. In higher dimensions we expect the use of walk on rectangles [9] 
techniques will facilitate the incorporation of complex geometries into our method, while still allowing the use of ba- 
sic Cartesian-grid meshes within each protective domain (hence avoiding the need for unstructured or embedded 
boundary meshes as commonly used for RDME models.) 

We have further demonstrated how our method can be utilized in practice. In particular, the many examples we 
investigated demonstrate that drift can have a significant effect on reaction locations, reaction time distributions, and 
number of reactions that occur in a system. To demonstrate how our approach might be used for more complicated 
systems, we investigated the process of diffusion-excursion of proteins interacting with a biopolymer (taking into ac- 
count geometric effects due to polymer shape). Such processes are thought to occur in the interactions between 
regulatory proteins and DNA. We considered a simplified model where the biopolymer has basic knotted configura- 
tions. Our results demonstrate that the combined effects of the biopolymer-protein interactions, biopolymer geometry, 
and drift-diffusion play a significant role in the chemical kinetics. The results also illustrate that capturing important 
features of biological systems may require models that account for drift-diffusion effects and spatial heterogeneities. 
We expect the FPKMC methods we have introduced to be useful in performing more realistic simulations of the 
chemical kinetics of biological systems. 
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A Derivation of Non-Uniform Jump Rates 

To derive the non-uniform rates in Eqn. (5), we follow the approach of Wang, Pesl<in, and Elston in [52] with the 
modification that the mesh width is no longer uniform. Let xi < xq < X2 or X2 < xq < xi be the locations of 
mesh points with non-uniform spacing hi = \xo — x\\ and = \xo - X2\, as illustrated in Fig. 3 { top row). We will 
denote the jump rates from xq to x^ by aoj, j = 1, 2. Let p{x, t) be the solution of the Fokl<er-Planck equation (3), 
which gives the probability density of being at location x at time t. Define Pi{t) to be the probability of being at the 
mesh point Xi at time t in the discrete model. Let p^^{x) and p^^ denote the equilibrium values of p{x, t) and Pi{t), 
respectively. We consider the point xi to represent the interval {xi - ^,xi + ^) \n the sense that 

Pi{t) ^ / p{x,t)dx p{xi,t)hi. 

Similarly, xq represents {xq- ^,xq + ^) and X2 represents {x2 - ^ ,X2 + 50 

Po (t) ^ p{xo,t) and p2 (t) ^ p{x2, t) /12 • 

We will look for jump rates aoj, from xo to Xj, by assuming that ajo, the jump rate in the opposite direction, is given 
by (4) with h = hj and imposing a discrete detailed balanced condition (zero net flux at equilibrium): 

aoj P°''(a:o) ^^ = ajoP°'^{xj)hj 

^ 2 D V{xo)-V{x,) p^%x,) 

"•^ /ii -I- /i2 hj exp[y(a;o) - V{xj)] - 1 p«^i(a;o) 

Since p^'^{x) is given by the Boltzmann distribution p'^'^{x) oc exp[-y(a;)], 

p-^x,) ^ exp[-y(x,)] ^ 1 

p«i(a;o) exp[-y(a;o)] exp[y(a;,) - y(a;o)] ' 

Thus, we obtain 

^ 2D V{x,) V{x,) 

°' hj{hi + h2)e^p[V{xj)-V{xo)]-l •' ' ■ 



B Analytic and Numerical Solutions for the Two-Molecule Annihilation Re- 
action A + B ^ 

In this Appendix, we discuss the numerical solution of the Fokker-Planck equation and the analytic solution of the 
diffusion equation to which we compared the FPKMC simulation results in Section 5.1 . The 1 D reaction-drift-diffusion 
system of two molecules undergoing the reaction A + B ^ can be described by Eqn. (6) on the 2D domain in 
the left panel of Figure 17. The solutions of the Fokker-Planck equation on the two disjoint triangular components of 
this domain are independent of each other. By symmetry, solving the Fokker-Planck equation on the domain in the 
left panel of Figure 17 can be reduced to solving the same equation on the single triangular domain in the center 
panel. We have written a PDE solver to solve the Fokker-Planck equation on the triangular domain. The PDF solver 
uses the rates (4) from the WPE discretization [52] of the Fokker-Planck equation. This discretization is second-order 
accurate for smooth potentials and first-order accurate for discontinuous potentials. The mesh in the PDE solver is 
uniform, and is chosen to be cell-centered at the Neumann boundaries and edge-centered at the Dirichlet boundary. 
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l = L-rR 1 = y/2{L - rn) 




Location of Molecule A Location of Molecule A x 

zero Dirichlet (absorbing, reactive) boundary _ _ _ zero Neumann (reflecting) boundary 

Figure 17: Left: 2D domain equivalent to 1 D simulation domain of length L in which the two molecules are located. 
The zero Dirichlet boundaries on the diagonals correspond to the reaction occurring between the two molecules when 
they are one reaction radius apart. The zero Neumann boundaries on the outer edges correspond to the reflecting 
endpoints of the 1 D simulation domain. Center: PDE solver domain. Right: Eigenfunction expansion domain. 



The solution to the diffusion equation on the triangular domain in Figure 17 (center) can be recovered by solving 
on a square domain formed by adjoining four copies of the triangular domain at their Neumann edges, as shown in 
the right panel of Figure 17. On the square domain, an eigenfunction expansion for the solution can be determined 
analytically. This provides a check for the FPKMC simulations in VneV — O case, as well as a check for the PDE 
solver in that case. The solution p{x,y,t) to the diffusion equation with constant initial condition po and diffusion 
coefficient D on a square domain with sides of length / is given by [41]: 



Wpo 



^ 2n + l 

n=Q 

oo 

E 



(2n + 1)ttx 
sm I ; I exp 



.m=0 



2m + 1 



(2to + l)7rw 
sm I I exp 



The survival probability S{t), the probability that the two particles have not yet reacted by time t, is 



S{t) 



I nl 



p{x, y, t) dx dy 



10 JQ 

The mean reaction time is 



E 

n=0 



1 



2n + 1 



exp 



/•oo /'OO fyA 740000^ 



-TT'^{2n + lfDt 



2m + 1 y (2n + 1)2 + (2m + 1)^ 



Evaluating E [T] at the parameter values corresponding to the two-molecule FPKMC simulations gives a mean 
reaction time of 0.06483188131 1 seconds. 

For each of the potential functions, we ran the PDE solver using the Crank-Nicolson method in time, with spatial 
step sizes ranging from Ax = tr = 0.02 to Aa; = rR/16 and time steps At ~ Ax/16. For F = 0, we can 
check the results of the PDE solver against the analytic solution; in this case, the numerical mean reaction times 
determined using the Crank-Nicolson method converge at approximately second-order (m = 2.0025) to the analytic 
mean reaction time. For V ^ there is no analytic solution to which the numerical solutions can be compared; 
however, the decrease in the pairwise differences between the mean reaction times determined from the numerical 
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Table 3: Improved errors in numerical survival probabilities. 





Crank-Nicolson 


TGA 




Aa; = rn/ie , Ai = Aa;/16 


Aa; = tr/IS , M = (Aa;)^ 


Error 


4.1117e-6 


4.79096-8 


Error 


8.25138-5 


1.61716-6 


L°° Error 


5.8634e-3 


2.2463e-4 



PDE solutions at successive step sizes indicates convergence (to = 1.9968 for V^os and to = 0.9844 for l^tep)- In 
the FPKMC convergence results in Section 5.1 , the empirical mean reaction times are compared to the analytic mean 
reaction time in the case V = 0, and to the numerical mean reaction times determined using the Crank-Nicolson 
PDE solver with the finest spatial step size, Aa; = tr/IG, in the cases of ycos and Kitep- 

The survival probabilities calculated using the Crank-Nicolson method show small numerical oscillations during 
the first few time steps. This is due to the incompatibility of the initial condition with the zero Dirichlet boundary 
condition. In order to numerically resolve the survival probabilities more accurately at short times, we re-ran the PDE 
solver using the Twizell-Gumel-Arigu (TGA) method [50] from t = to t = 0.07 seconds. The TGA method is a 
second order, Lq stable time discretization. We also used a finer time step, M = (Aa;)^ where Aa; = tr/IG, to 
further improve the accuracy at short times when the survival probabilities change most rapidly. We then determined 
the numerical survival probabilities using the results from the TGA method for f = to t = 0.07 seconds and using 
the results from Crank-Nicolson method for t > 0.07 seconds. Att = 0.07, the absolute difference in survival 
probabilities between the two methods is less than 10~^ for V^cos and less than 10"'^ for Vstep- 

To check that using the TGA method with a finer time step improved the accuracy of the numerical survival 
probabilities over the Crank-Nicolson method, we compared the numerical survival probabilities for y = to the 
analytic survival probability on the interval from t = OXot = 0.07. Table 3 shows the absolute errors between the 
numerical survival probability calculated using each method and the analytic survival probability. 
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